# install package librarian if neededif (!("librarian"%in%rownames(installed.packages()))) {install.packages("librarian")}# load required packageslibrarian::shelf( tidyverse, usmap, fs, ggpubr, ggdist, ggrepel, faux, lme4, lmerTest, ggeffects, binom, tictoc, ggthemes, sessioninfo, knitr, kableExtra)# Source required functionsmyFunctions <-c("FUNStormEventsData_filterData")for (f in myFunctions) {source(paste0("../functions/", f, ".R"))}# Preperations to show states boundariespoly_states <-plot_usmap(regions ="states")# Read in data_details_fipsfileName <-"data_details_fips.RDS"pathName <-"../data/stormData"filePath <-dir_ls(path = pathName, regexp =paste0(fileName, "$")) %>%last()data_details_fips <-readRDS(filePath)
1 Carbon Emission Task: Trials
Table 1 shows all trials in the new variant of the Carbon Emission Task (CET). Each row corresponds to one trial. Carbon Diff. and Bonus Diff. index the unique combinations of relative differences in carbon and bonus consequences in a trial. Carbon Level Rand. and Bonus Level Rand. correspond to the random carbon and bonus level used in the construction of a trial. The average carbon level over all trials is 19.85 lbs CO2, to which we add a random deviation with mean = 0 and SD = 0.9925 lbs CO2. The average bonus level over all trials is 60 cents, to which we add a random deviation with mean = 0 and SD = 3 cents USD. Carbon Option Pro-Self and Bonus Option Pro-Self represent the actual information available to participants related to the option that will lead to the greater bonus payment for participants. Carbon Option Pro-Climate and Bonus Option Pro-Self represent the same information but for the option leading to the greater benefit for the climate (lower carbon emissions).
Show the code
# Read in the trials dataCETTrials <-read_csv("../data/CET_trials.csv")# Select columns to display and round to two digitsCETTrials <- CETTrials %>%select( carbon_prcnt, bonus_prcnt, carbon_level_rand, bonus_level_rand, carbon_self, carbon_env, bonus_self, bonus_env ) %>%mutate(across(-ends_with("_prcnt"), ~round(.x, 2)))# Define column namesCETTrials_colnames <-c("Carbon Diff (%)","Bonus Diff (%)","Carbon Level Rand (lbs.CO2)","Bonus Level Rand (cent USD)","Carbon Option Pro-Self (lbs.CO2)","Carbon Option Pro-Climate (lbs.CO2)","Bonus Option Pro-Self (cent USD)","Bonus Option Pro-Climate (cent USD)")# Dispaly table with custom width for certain columnskable(CETTrials, col.names = CETTrials_colnames) %>%column_spec(column =1:2, width ="1cm") %>%column_spec(column =3:8, width ="2cm")
Table 1: Trials in the new variant of the CET
Carbon Diff (%)
Bonus Diff (%)
Carbon Level Rand (lbs.CO2)
Bonus Level Rand (cent USD)
Carbon Option Pro-Self (lbs.CO2)
Carbon Option Pro-Climate (lbs.CO2)
Bonus Option Pro-Self (cent USD)
Bonus Option Pro-Climate (cent USD)
10
10
20.54
60.33
21.56
19.51
63
57
10
15
19.80
63.84
20.79
18.81
69
59
10
20
19.89
57.64
20.89
18.90
63
52
10
50
20.78
58.24
21.82
19.74
73
44
10
100
20.44
62.12
21.46
19.42
93
31
15
10
20.14
59.66
21.66
18.63
63
57
15
15
19.37
56.33
20.83
17.92
61
52
15
20
18.69
63.63
20.09
17.29
70
57
15
50
19.47
63.24
20.93
18.01
79
47
15
100
18.37
64.15
19.75
16.99
96
32
20
10
19.27
60.23
21.20
17.34
63
57
20
15
18.55
63.58
20.40
16.69
68
59
20
20
19.31
63.65
21.24
17.38
70
57
20
50
19.56
60.05
21.52
17.60
75
45
20
100
20.38
63.13
22.42
18.34
95
32
50
10
19.30
57.58
24.13
14.48
60
55
50
15
20.82
55.45
26.02
15.61
60
51
50
20
20.46
59.72
25.58
15.35
66
54
50
50
20.89
57.35
26.12
15.67
72
43
50
100
21.21
60.94
26.52
15.91
91
30
100
10
19.89
60.80
29.83
9.94
64
58
100
15
20.54
61.23
30.81
10.27
66
57
100
20
19.54
65.30
29.30
9.77
72
59
100
50
20.50
60.10
30.75
10.25
75
45
100
100
19.04
63.07
28.57
9.52
95
32
2 Extreme Weather Data
2.1 Purpose & Rationale
As outlined in the Registered Report, we will assess the number of extreme weather episodes recorded in each participant’s county of residence within the 30 days prior to study completion. Regarding the time window during which we plan to conduct the study, we aim for maximizing the likelihood of capturing suitable variability in the exposure to extreme weather episodes with notable geographic variability. To this end, we analyzed records of extreme weather episodes over the last ten years.
2.2 Filter Data
We filter the storm events data for the specific years, months, and extreme weather event types we are interested in. We filter for all years from 2014 to 2023 (as data are not complete for the year 2024 yet), we highlight the month of July, and we focus on those types of extreme weather events that are predicted to increase in frequency and severity due to climate change : Excessive Heat, Drought, Wildfire, Flash Flood, Coastal Flood, Strong Wind, Hail, and Tornado (IPCC 2023).
Our analysis spanning the last ten years evidenced high numbers of occurrences of extreme weather episodes in every single year (Figure 1) as well as considerable geographical variability (Figure 2). Analyzing the seasonal distribution of extreme weather episodes, Figure 1 shows that July consistently shows a high number of extreme weather episodes over the last ten years, with an average occurrence of 1,777 episodes in each year (SD = 450). Note that numbers represent counts of distinct episode-county combinations to represent how many counties were affected. Table 2 shows the average number of extreme weather events by event type in July over the last ten years, demonstrating that July typically witnesses a multitude of different extreme weather event types. Additionally, Figure 2 indicates that within the month of July, these extreme weather episodes also display a high geographical variability.
Figure 1: Histograms showing the number of extreme weather episodes by month from 2014 to 2023. The dashed horizontal line indicates the mean number of extreme weather episodes in each year. The thick-bordered bar marks the month with the most extreme weather events each year. The orange bar represents July. July had the most extreme weather events in 4 out of 10 years, and in another 4 years, it was right before or after the peak month. Only episodes that included at least one of the following event types were considered: excessive heat, drought, wildfire, flash flood, coastal flood, strong wind, hail, tornado.
Table 2: Mean Number of Extreme Weather Events in July over the Years 2014 to 2023 by Event Type.
Event Type
Mean Number of Events
SD
Coastal Flood
5.43
3.00
Drought
272.70
254.00
Excessive Heat
429.50
347.00
Flash Flood
778.80
241.00
Hail
1249.70
383.00
Strong Wind
8.60
6.00
Tornado
107.60
20.00
Wildfire
106.90
34.00
Figure 2: Maps displaying the geographical distribution of the occurrence of at least one extreme weather episode in July over the years 2014 to 2023. Only episodes that included at least one of the following event types were considered: excessive heat, drought, wildfire, flash flood, coastal flood, strong wind, hail, tornado.
Show the code
dataForPlot <- out$dataForUsPlot %>%mutate(nEpisodes_withNA =ifelse(nEpisodes ==0, NA_integer_, nEpisodes))p.map_cont <-plot_usmap(data = dataForPlot,values ="nEpisodes_withNA",regions ="counties",exclude =c("AK", "HI"),color ="black",linewidth =0.1 ) +geom_sf(data = poly_states[[1]] %>%filter(!(abbr %in%c("AK", "HI"))),color ="black",fill =NA,linewidth = .3 ) +scale_fill_binned(name ="Number of Episodes",n.breaks =10,type ="viridis",na.value ="white" ) +labs(title ="Extreme Weather Episodes in July over the Years 2014 to 2023" ) +theme_bw() +theme(text =element_text(size =15),legend.position ="bottom",plot.title =element_text(hjust = .5),panel.grid =element_blank(),axis.ticks =element_blank(),axis.text =element_blank() ) +facet_wrap(~year, ncol =5)jpeg(file ="../images/mapGeographicalDistribution_cont.jpeg",width =14, height =7.5, units ="in", res =600)print(p.map_cont)invisible(dev.off())p.hist_count <- out$dataForUsPlot %>%group_by(year, nEpisodes) %>%summarise(count =n(),prcnt = count /n_distinct(out$dataForUsPlot$fips) ) %>%ggplot(aes(x = nEpisodes, y = prcnt)) +geom_bar(stat ="identity", color ="black", fill ="darkgrey") +scale_y_continuous(labels = scales::label_percent()) +labs(x ="Number of Episodes",y ="Proportion of Counties" ) +theme_bw() +labs(title ="Extreme Weather Episodes in July over the Years 2014 to 2023" ) +theme(text =element_text(size =15),legend.position ="bottom",plot.title =element_text(hjust = .5) ) +facet_wrap(~year, ncol =5)jpeg(file ="../images/frequencyDistribution_cont.jpeg",width =14, height =7.5, units ="in", res =600)print(p.hist_count)invisible(dev.off())# Calcualte some proportions for display in textprops2023 <- out$dataForUsPlot %>%filter(year ==2023) %>%count(episodes_bin) %>%mutate(freq = n/sum(n),freq_prcnt =paste0(format(round(freq*100, 2), nsmall =2), "%") )
While Figure 2 visualizes the occurrence of at least one extreme weather episode in July for each county and year (binary variable), Figure 4 displays the actual number of such episodes (continuous). The vast majority of counties were exposed to few episodes, indicating that most of the variability is due to whether an extreme weather episode occurred at all or not. This is further supported by Figure 3 showing histograms for the number of extreme weather episodes in July over the past ten years. Most counties reported either zero or one extreme weather episode in July, and the ratio of counties experiencing no episodes to counties experiencing at least one episode seems to gradually approach 1:1. In July 2023, for instance, this ratio reached 1.02, with 50.43% of counties being exposed to zero and 49.57% of counties being exposed to at least one extreme weather episode.
Figure 3: Maps displaying the geographical distribution of the raw number of extreme weather episodes in July over the years 2014 to 2023. The color palette indicates numbers greater than zero, and white represent a count of zero episodes.
Figure 4: Histograms displaying the distribution of the raw number of extreme weather episodes in July over the years 2014 to 2023. For each number of episodes on the x-axis, the y-axis shows the proportion of counties that recorded this number of episodes.
Finally, as reported in the analysis plan and the design table, we plan to run a set of additional analyses regarding hypotheses H2 and H3, in which we will test the sensitivity of results to the time period prior to study completion used to assess extreme weather exposure. Regarding H2, we will estimate the two-way interaction effect of political affiliation and extreme weather exposure on ΔDuration for different time periods from 30 days to 360 days in increments of 30 days. Similarly for H3, we will estimate the three-way interaction effect of political affiliation, extreme weather exposure, and attribution of extreme weather events to climate change on ΔDuration for the same time periods. We will visualize results of these additional analyses by plotting the two-way (or three-way) interaction regression coefficients as points surrounded by their 95%-CI on the y-axis and the 12 time periods on the x-axis, as displayed in Figure 5 with simulated data. Based on previous research (Konisky, Hughes, and Kaylor 2016), we expect that the estimated effects will decay as the number of days prior to study completion used to assess the occurrence of extreme weather episodes increases.
Figure 5: Simulated regression coefficients of interaction effects for different number of days prior to study completion used to assess the occurrence of extreme weather episodes. Point estimates are surrounded by their 95% confidence intervals. The dashed line represents the absence of an interaction effect (regression coefficient of zero). Significance of regression coefficients is color-coded, with black indicating regression coefficients significantly different from zero, and grey indicating no significant difference from zero.
2.4 Conclusion
Our analyses indicate that July consistently shows a high number of extreme whether episodes with notable geographic variability (Figure 1 and Figure 2). Therefore, to maximize the likelihood of capturing suitable variability in exposure to extreme weather episodes, we plan to conduct our study at the beginning of August, ensuring that the 30-day period prior to study completion falls within July. Moreover, the main source of variability in exposure to extreme weather episodes in July is due to whether at least one episode occurred or not (Figure 4 and Figure 4). Thus, our main analyses will focus on whether a participant was exposed to at least one extreme weather episode in the 30 days prior to study completion, treated as a binary variable. In additional analyses, we will test the sensitivity of our results to different time periods used to assess extreme weather exposure prior to study completion.
3 Power Analyses
3.1 Purpose & Rationale
The primary goal of the present analyses is to determine the number of participants required to achieve 95% statistical power to detect a Smallest Effect Size Of Interest (SESOI) for a main effect of political affiliation (Democrat vs. Republican) on attentional information search behavior as assessed by ΔDuration. That is, we primarily aim for conducting a power analysis for sample-size determination, also called a priori power analysis (Giner-Sorolla et al. 2024). To this end, we use power simulations with parameters informed by previous studies to assess the statistical power for different sample sizes. This will allow us to decide on a sample size that will provide high statistical power to detect a true and theoretically relevant main effect of political affiliation.
Given this sample size, our secondary goal is to then assess the statistical power to detect different effect sizes for (1) the two-way interaction between political affiliation and extreme weather exposure and (2) the three-way interaction between political affiliation, extreme weather exposure, and the subjective attribution of extreme weather events to climate change. That is, we conduct power analyses for assessing effect-size sensitivity for these two- and three-way interactions (Giner-Sorolla et al. 2024).
The present report is organized as follows:
We describe all important variables of the planned study and how we will assess them.
We derive a SESOI for the main effect of political affiliation on ΔDuration.
We conduct power simulations for sample-size determination analysis based on this political affiliation main effect SESOI.
We conduct power simulations for effect-size sensitivity analyses for a two-way interaction effect between political affiliation and extreme weather exposure.
We conduct power simulations for effect-size sensitivity analyses for a three-way interaction effect between political affiliation, extreme weather exposure, and subjective attribution of extreme weather events to climate change.
3.2 Important Variables
3.2.1 ΔDuration
Participants will complete 25 trials in a new variant of the Carbon Emission Task (CET, Berger and Wyss 2021) optimized for online process-tracing using MouselabWEB (Willemsen and Johnson 2019). An example trial of the task is displayed and explained in Figure 6.
Figure 6: An example trial of the new variant of the CET. (A) Participants are presented with two options which are associated with different bonus payment and carbon emission consequences. (B) Participants inspect the exact consequences regarding each attribute of each option by hovering their mouse over the respective boxes. Moving the mouse outside of a box occludes the information again. Note that whether the bonus or carbon attributes are displayed in the top row is randomized across participants but held constant within participants. (C) By visiting all boxes, participants can acquire all available information in a trial (all boxes opened simultaneously for demonstration purposes only). Note that whether the option maximizing the bonus payment for participants is presented in the left or right row is randomized within participants.
In all analyses, the criterion (dependent variable) will be ΔDuration (deltaDuration), which is calculated in each trial as:
with \(t\) representing the summed up dwell time on \(Carbon/Bonus\) boxes in Option \(A/B\). ΔDuration varies between -1 and 1 with:
a value of -1 indicating that the entire dwell time was spent on gathering Bonus information
a value of 0 indicating that the dwell time was equally split between gathering Bonus and Carbon information
a value of 1 indication that the entire dwell time was spent on gathering Carbon information
3.2.2 Political Affiliation
We will assess political affiliation using the following question:
Generally speaking, do you think of yourself as a Democrat, Republican or Independent?
Participants will answer on the following 7-point Likert scale:
[1] Strong Republican, [2] Not strong Republican, [3] Independent, close to Republican, [4] Independent, [5] Independent, close to Democrat, [6] Not strong Democrat, [7] Strong Democrat
Additionally, if participants self-identify as [4] Independent, they will be asked:
You said that you think of yourself as an Independent politically. If you had to identify with one party of the two parties, which one would you choose?
Participants will answer on the following scale:
[1] Republican Party, [2] Democratic Party
Based on these questions, we classify participants as either Republican or Democrat as represented by the variable polAff that can take on the following values:
[rep] Republican Party, [dem] Democratic Party
3.2.3 Extreme Weather Exposure
For each participant, we will assess the number of extreme weather episodes that occurred in the participants’ county of residence in the 30 days prior to study completion. We then create a variable ewe which equals to TRUE if at least one extreme weather episode occurred in the specified time interval, and FALSE otherwise.
3.2.4 Subjective Attribution
We will assess the degree to which participants attribute extreme weather events to climate change (see Ogunbode et al. 2019). Participants will rate their agreement to the following three questions:
Extreme weather events are caused in part by climate change.
Extreme weather events are a sign that the impacts of climate change are happening now.
Extreme weather events show us what we can expect from climate change in the future.
Participants will answer on the following 5-point Likert scale:
Subjective attribution of EWE to climate change will be operationalised as the mean agreement to these three statements. Note that Ogunbode et al. (2019), who used the same questions, response options, and aggregation, report a mean of 3.67 and a SD of 0.85 for this variable.
3.3 SESOI
As argued in the Registered Report, we hypothesize that (H1) compared to Republicans, Democrats prioritize searching for and attending to carbon over bonus information during decision-making in the CET. That is, Democrats display higher ΔDuration values compared to Republicans.
While H1 describes the direction of the expected effect, it does not specify the expected magnitude of the effect. This later point is clarified by asking the question what would be the smallest effect size that researchers would still consider theoretically relevant, i.e., what is the smallest effect size of interest (SESOI)? We argue for such a SESOI on theoretical grounds and based on previous MouselabWEB research.
In mouselabWEB studies, it is standard practice to filter out any information acquisition lasting less than 200 msec because such very short (spurious) acquisitions are very unlikely to be consciously processed (Willemsen and Johnson 2019). Therefore, we derive the SESOI based on the consideration that for an effect to be meaningful, the increase in ΔDuration from a Republican to a Democrat should be due to an increase of the time spent gathering carbon (relative to bonus) information of at least 200 msec. We need to translate these considerations into the metric of ΔDuration.
Suppose that Republicans on average spend \(t_{Carbon}\) msec on gathering carbon information and \(t_{Bonus}\) msec on gathering bonus information in each trial, with a total time of gathering any information \(t_{Total} = t_{Carbon} + t_{Bonus}\). Thus, for Republicans, this results in:
Suppose that Democrats have the same \(t_{Total}\) as Republicans, but they spend 200 msec more on gathering carbon information and, correspondingly, 200 msec less on gathering bonus information. Thus, for Democrats, this results in:
As this definition shows, the SESOI depends on our expectation regarding \(t_{Total}\). We form these expectations based on previous research. In our lab, we recently conducted another MouselabWEB study whose setup was very similar to the one described in the Registered Report (Studler et al., in preparation). In short, student participants completed a behavioral task in which they searched for and attended to information presented in a 2-by-2 decision matrix adapted from Reeck, Wall, and Johnson (2017). The average time participants spent on acquiring information in each trial (\(t_{Total}\)) was 3.3 seconds. In the planned study we will assess a representative US sample. That is, the average age of our sample will be higher than in a typical student sample. As processing speed is known to decline with age (Salthouse 2000), we assume a slightly higher total information acquisition time in our sample of \(t_{Total} = 3500\ msec\). Therefore, we define our SESOI as:
To account for uncertainty in this estimation, we also provide sample-size determination analysis results based on \(t_{Total} = 4000\ msec\) and \(t_{Total} = 3000\ msec\). This results in a high and low estimate of SESOI:
As argued in the Registered Report, we hypothesize:
(H1)Compared to Republicans, Democrats prioritize searching for and attending to carbon over bonus information during decision-making in the Carbon Emission Task. That is, Democrats display higher ΔDuration values compared to Republicans.
Show the code
# Create a data frame with predicted true effects# Smallest Effect Size Of Interest (SESOI)SESOI <-0.4/3.5# Betasbeta_p <- SESOI# Predicted "true" effectspolAff_trueEffects <-expand_grid(polAff =factor(c("rep", "dem"), levels =c("rep", "dem"))) %>%add_contrast("polAff", contrast ="anova", colnames ="X_p") %>%mutate(trueDeltaDuration =0+# intercept X_p * beta_p # main effect polAff )
3.4.1 Data Simulation Function
We first define a function that simulates data for the main effect of political affiliation on ΔDuration: FUN_sim. The function will simulate data according to the following model, expressed in lme4-lingo:
deltaDuration ~ polAff + (1|subj) + (1|trial)
The function FUN_sim takes, among others, the following important arguments:
n_subj: Number of subjects. Chaning this parameter allows us to assess statistical power for different sample sizes.
beta_0: Fixed intercept (grand mean). We assume that the average participant (irrespective of political affiliation or any other predictor variable) spends about the same time on searching for and attending to carbon as to bonus information in the CET. Therefore, we set beta_0 to zero. The effects of predictor variables will be modeled as deviations from this grand mean.
beta_p: Fixed effect of political affiliation. This value represents the average difference in ΔDuration between Democrats and Republicans (Democrats - Republicans). As discussed above, we set this value to 0.1143 by default, and we assess the effect of changing this variable on statistical power.
subj_0: By-subject random intercept SD. We simulate that a participants’ deviations from the grand mean for ΔDuration follows a normal distribution with a mean of 0 and a standard deviation of subj_0 = 0.29. We base our default value on a study by Reeck, Wall, and Johnson (2017). They investigated whether variability in information search behavior is driven predominantly by differences in the features of a choice (i.e., in our case: the relative differences between carbon and bonus outcomes in options A and B) or by individual differences. To this end, they predicted information acquisition using a intercept-only model that included random intercepts for subjects and items. They estimated the random intercept of subjects to be 0.29.
trial_0: By-trial random intercept SD. We simulate that a items’ deviations from the grand mean for ΔDuration follows a normal distribution with a mean of 0 and a standard deviation of trial_0 = 0.04. Again, we base our default value on the study by Reeck, Wall, and Johnson (2017), who estimated a random intercept for trials of 0.04.
sigma: Trial-level noise (error) SD. We model the error SD to be of the same size as the sum of the random intercept SDs = 0.29 + 0.04 = 0.33. We also report simulation results for an error SD that is twice the size of the random intercept SDs, i.e., 0.66.
The function FUN_sim is defined below (code visible in html output only):
Show the code
# define data simulation functionFUN_sim <-function(n_subj =1000, # number of subjectsn_subj_prop =c(.5, .5), # proportion of republican and democrat subjectsn_trial =25, # number of trialsbeta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # effect of political affiliation on deltaDurationsubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible deltaDuration values be truncuated?setSeed =NULL# seed number to achieve reproducible results. Set to NULL for simulations!) {# set seed to achieve reproducible results for demonstration purposesset.seed(setSeed)# simulate data for dwell time on carbon information dataSim <-# add random factor subjectadd_random(subj = n_subj) %>%# add random factor trialadd_random(trial = n_trial) %>%# add between-subject factor political affiliation (with anova contrast)add_between("subj", polAff =c("rep", "dem"), .prob = n_subj_prop*n_subj, .shuffle =FALSE) %>%add_contrast("polAff", colnames ="X_p", contrast ="anova") %>%# add by-subject random interceptadd_ranef("subj", S_0 = subj_0) %>%# add by-trial random interceptadd_ranef("trial", T_0 = trial_0) %>%# add error termadd_ranef(e_st = sigma) %>%# add response valuesmutate(# add together fixed and random effects for each effectB_0 = beta_0 + S_0 + T_0,B_p = beta_p,# calculate dv by adding each effect term multiplied by the relevant# effect-coded factors and adding the error termdeltaDuration = B_0 + (B_p * X_p) + e_st )# unset seedset.seed(NULL)# truncuate impossible deltaDurationsif(truncNums) { dataSim <- dataSim %>%mutate(deltaDuration =if_else(deltaDuration <-1, -1,if_else(deltaDuration >1, 1, deltaDuration))) }# run a linear mixed effects model and check summary mod <-lmer( deltaDuration ~ polAff + (1| subj) + (1| trial),data = dataSim, ) mod.sum <-summary(mod)# get results in tidy format mod.broom <- broom.mixed::tidy(mod)return(list(dataSim = dataSim,modelLmer = mod,modelResults = mod.broom ))}
We call the function once and extract the results of this single simulation (code visible in html output only):
Show the code
# Call functionout <-FUN_sim(n_subj =1000, # number of subjectsn_subj_prop =c(.5, .5), # proportion of republican and democrat subjectsn_trial =25, # number of trialsbeta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # effect of political affiliation on deltaDurationsubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible deltaDuration values be truncuated?setSeed =1234# seed number to achieve reproducible results. Set to NULL for simulations!)# Get results tableresultsTable <- out$modelResults %>%select(-c(std.error, statistic, df)) %>%mutate(across(where(is_double), ~round(.x, 4))) %>% knitr::kable()formulaUsedForFit <-paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse =" ")# Create plotp.demo.mainEffect <- out$dataSim %>%ggplot(aes(x = polAff, y = deltaDuration, color = polAff)) +geom_hline(yintercept =0) +geom_violin(alpha =0.3) +geom_point(data = polAff_trueEffects,mapping =aes(x = polAff, y = trueDeltaDuration),shape ="circle open",size =3.5,stroke =2,color ="black" ) +stat_summary(fun = mean,fun.min = \(x){mean(x) -sd(x)},fun.max = \(x){mean(x) +sd(x)},position =position_dodge(width = .9) ) + ggrepel::geom_label_repel(data = polAff_trueEffects,mapping =aes(x = polAff, y = trueDeltaDuration, label =round(trueDeltaDuration, 4)),color ="black",box.padding =1 ) +scale_color_manual(values =c("red", "dodgerblue")) +scale_y_continuous(breaks =seq(-1, 1, .2)) +labs(title ="Demo Output of One Simulation for Main Effect") +theme_bw() +theme(legend.position ="none")jpeg(file ="../images/pa_demoMainEffect.jpeg",width =9.1, height =6.5, units ="in", res =600)print(p.demo.mainEffect)invisible(dev.off())
Figure 7 visualizes the results of this single simulation and Table 3 summarises the statistical results of fitting the actual model used in data generation to the simulated data.
Figure 7: Visual representation of results of one simulation created using FUN_sim. Violin plots display the full distribution of the data. Points and surrounding lines indicate the mean ± 1 SD. The black horizontal line displays the true sample mean and the black open circles indicate the true means for each cell.
Show the code
resultsTable
Table 3: Statistical results of one simulation created using FUN_sim. Data was fit using deltaDuration ~ polAff + (1 | subj) + (1 | trial).
effect
group
term
estimate
p.value
fixed
NA
(Intercept)
-0.0188
0.1105
fixed
NA
polAff.dem-rep
0.0896
0.0000
ran_pars
subj
sd__(Intercept)
0.2841
NA
ran_pars
trial
sd__(Intercept)
0.0363
NA
ran_pars
Residual
sd__Observation
0.3223
NA
3.4.2 Power Simulation
We now simulate multiple random samples drawn from the same synthetic population with a known true effect of political affiliation. For each random sample, we fit our statistical model to the data. The statistical power to detect a true effect of political affiliation is calculated as the proportion of significant effects out of the total number of simulations. We aim for a statistical power of 95%.
In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run (code only shown in the html version).
Show the code
FUN_sim_pwr <-function(sim, ...){ out <-FUN_sim(...) modelResults <- out$modelResults %>%mutate(sim = sim) %>%relocate(sim)return(modelResults)}# How many simulations should be run?n_sims <-1000# What are the breaks for number of subjects we would like to calculate power for?breaks_subj <-seq(200, 1000, 200)# What are the breaks for SESOI?breaks_sesoi <-c(0.4/3, 0.4/3.5, 0.4/4)# What are the breaks for different errer SDs?breaks_sigma <-c(1:2)*(.29+.04)res_mainEffect <-tibble()for (s inseq_along(breaks_sigma)) { res_sesoi <-tibble()for (sesoi inseq_along(breaks_sesoi)) { res_nSubj <-tibble()for (nSubj inseq_along(breaks_subj)) {# Give feedback regarding which model is simulatedcat(paste0("Simulation:\n"," sigma = ", round(breaks_sigma[s], 4), "\n"," sesoi = ", round(breaks_sesoi[sesoi], 4), "\n"," nSubject = ", breaks_subj[nSubj], "\n" ))# Start timercat(paste0("Start date time: ", lubridate::now(), "\n"))tic()# Loop over simulations pwr <-map_df(1:n_sims, FUN_sim_pwr,n_subj = breaks_subj[nSubj],beta_p = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Stop timer and calculate elapsed time elapsed_time <-toc(quiet =TRUE) elapsed_seconds <- elapsed_time$toc - elapsed_time$tic elapsed_minutes <- elapsed_seconds /60cat(paste0("End date time: ", lubridate::now(), "\n"))cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")# Add number of subjects to pwr pwr <- pwr %>%mutate(nSubjects = breaks_subj[nSubj],sesoi = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Add results to the results table res_nSubj <- res_nSubj %>%rbind(pwr) }# Add results to the results table res_sesoi <- res_sesoi %>%rbind(res_nSubj) }# Add results to the results table res_mainEffect <- res_mainEffect %>%rbind(res_sesoi)}res_mainEffect.summary <- res_mainEffect %>%filter(term =="polAff.dem-rep") %>%group_by(sigma, sesoi, nSubjects) %>%summarise(power =mean(p.value <0.05),ci.lower =binom.confint(power*n_sims, n_sims, methods ="exact")$lower,ci.upper =binom.confint(power*n_sims, n_sims, methods ="exact")$upper,.groups ='drop' ) %>%mutate(sigma_fact =factor(format(round(sigma, 4), nsmall =4)),sigma_level =match(sigma_fact, levels(sigma_fact)),sesoi_fact =factor(format(round(sesoi, 4), nsmall =4)),sesoi_level =match(sesoi_fact, levels(sesoi_fact)) )# Save results in a list objecttime <-format(Sys.time(), "%Y%m%d_%H%M")fileName <-paste0("res_mainEffect", "_", time, ".RDS")saveRDS(list(res_mainEffect = res_mainEffect,res_mainEffect.summary = res_mainEffect.summary ),file =file.path("../powerSimulationsOutput", fileName))
We retrieve pre-run results:
Show the code
# Load power simulation dataresList_mainEffect <-readRDS(file.path("../powerSimulationsOutput", "res_mainEffect_20240813_2141.RDS"))resList_mainEffect.summary <- resList_mainEffect$res_mainEffect.summary# Extract power values for some specific effect sizes at N = 1000chosenN <-600chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1000", "0.1143")powerValues <- resList_mainEffect.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN) %>%mutate(power_str =paste0(round(ci.lower*100, 2), "%")) %>%pull(power_str)# Extract number of simulationslabel_nSimulations <- resList_mainEffect$res_mainEffect$sim %>%n_distinct()
Figure 8 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.
Figure 8: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.
Figure 9 shows results of our sample-size determination analyses. We find that a sample size of 600 provides a statistical power of 95.4% (lower bound of 95%-CI) even under the most conservative assumptions (SESOI = 0.1000, Error SD = 0.6600).
Show the code
# Get power values of interestchosenN <-c(800, 1000)chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1000")powerValues <- resList_mainEffect.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN)# Get lower CI for power values for conservative sigma assumptionspowerValues_sigmaCons <- powerValues %>%filter(sigma_fact =="0.6600")# Interpolate the effect sizes at which we achieve 99% powereff_sigmaCons <-approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$nSubjects, xout =0.99)$y# Round results for display in texteff_sigmaCons_txt <-toString(round(eff_sigmaCons))
Moreover, we inspect the lower bounds of the 95%-CI power estimates in Figure 9. Specifically, we focus on the power simulation results for a conservative SESOI of 0.1000 and an error SD of 0.6600. Then, we use interpolation between sample sizes of 800 and 1000 to estimate the number of participants needed to achieve 99% statistical power (lower bounds of 95%-CI power estimates). Results reveal that with a sample size of N = 895, we would achieve 99% statistical power (lower bound of 95%-CI power estimate) to detect a true effect of at least 0.1000.
We optimize the study design to detect a true SESOI for political affiliation. However, we are also interested in two-way and three-way interaction effects, which are known to require greater sample sizes to achieve the same statistical power as for main effects. Moreover, greater sample sizes are more likely to accurately represent target populations with respect to variables like exposure to extreme weather events and subjective attribution of extreme weather events to climate change. Therefore, we opt for a sample size of N = 1000 for further effect-size sensitivity analyses regarding interaction effects.
Figure 9: Power curves for the main effect of polAff. Points represent statistical power surrounded by a 95%-CI based on 1000 simulations with α = 0.05.
3.5 Two-Way Interaction Effect
As argued in the Registered Report, we hypothesize:
(H2): Exposure to extreme weather events amplifies the polarization in attentional information search processes, increasing the difference in ΔDuration between Republicans and Democrats among individuals with (vs. without) such exposure. In other words, there is a positive two-way interaction effect of extreme weather exposure and political affiliation on ΔDuration.
For deriving a SESOI for the two-way interaction of interest, similar considerations apply as in the case of the SESOI for the main effect of interest. In this section, we make these considerations explicit. We start by noticing that the complete fixed two-way interaction polAff × ewe is modeled as:
Now, let’s calculate this effect for two individuals who differ in their levels of ewe. As ewe is a binary variable, we have two types of individuals: individuals with (ewe = 1) and without (ewe = 0) extreme weather exposure. For an individual without extreme weather exposure, the effect of polAff will be:
We next define a function that simulates data for the two-way interaction effect of political affiliation with extreme weather exposure on ΔDuration: FUN_sim_2wayInt. The function will simulate data according to the following model, expressed in lme4-lingo:
The function FUN_sim_2wayInt takes, among others, the following important arguments (in addition to the arguments discussed for FUN_sim):
beta_p: Fixed main effect of political affiliation. Compatible with H1, we keep this value at the SESOI of 0.1143. That is, we model that the average effect of political affiliation across all participants, irrespective of their extreme weather exposure, is 0.1143.
beta_e: Fixed main effect of extreme weather exposure. As we are interested in the moderating role of extreme weather exposure, we set this main effect to zero. That is, we assume that the effect of extreme weather exposure is highly dependent on participants’ political affiliation.
beta_p_e_inx: Fixed two-way interaction effect of political affiliation and extreme weather exposure. We set this initial value to the SESOI derived above, but we investigate how changing this effect size impacts statistical power, as we are conducting effect-size sensitivity analyses for interaction effects.
The function FUN_sim_2wayInt is defined below:
Show the code
# define data simulation functionFUN_sim_2wayInt <-function(n_subj =1000, # number of subjectsn_subj_prop_p =c(.5, .5), # proportion of republican and democrat subjectsn_subj_prop_e =c(.5, .5), # proportion of subjects without and with extreme weather exposuren_trial =25, # number of trialsbeta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # main effect of political affiliation (polAff)beta_e =0, # main effect of extreme weather exposure (ewe)beta_p_e_inx =0.4/3.5, # two-way interaction effect of polAff and ewesubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible numbers be truncuated?setSeed =NULL# seed number to achieve reproducible results. Set to NULL for simulations!) {# set seed to achieve reproducible results for demonstration purposesset.seed(setSeed)# simulate data for dwell time on carbon information dataSim <-# add random factor subjectadd_random(subj = n_subj) %>%# add random factor trialadd_random(trial = n_trial) %>%# add between-subject factor political affiliation (with anova contrast)add_between("subj", polAff =c("rep", "dem"), .prob = n_subj_prop_p*n_subj, .shuffle =TRUE) %>%add_contrast("polAff", colnames ="X_p", contrast ="anova") %>%# add between-subject factor extreme weather exposure (with anova contrast)add_between("subj", ewe =c(FALSE, TRUE), .prob = n_subj_prop_e*n_subj, .shuffle =TRUE) %>%add_contrast("ewe", colnames ="X_e", contrast ="anova") %>%# add by-subject random interceptadd_ranef("subj", S_0 = subj_0) %>%# add by-trial random interceptadd_ranef("trial", T_0 = trial_0) %>%# add error termadd_ranef(e_st = sigma) %>%# add response valuesmutate(# add together fixed and random effects for each effectB_0 = beta_0 + S_0 + T_0,B_p = beta_p,B_e = beta_e,B_p_e_inx = beta_p_e_inx,# calculate dv by adding each effect term multiplied by the relevant# effect-coded factors and adding the error termdeltaDuration = B_0 + e_st + (X_p * B_p) + (X_e * B_e) + (X_p * X_e * B_p_e_inx) )# truncuate impossible deltaDurationsif(truncNums) { dataSim <- dataSim %>%mutate(deltaDuration =if_else(deltaDuration <-1, -1,if_else(deltaDuration >1, 1, deltaDuration))) }# run a linear mixed effects model and check summary mod <-lmer( deltaDuration ~ polAff*ewe + (1| subj) + (1| trial),data = dataSim ) mod.sum <-summary(mod)# get results in tidy format mod.broom <- broom.mixed::tidy(mod)return(list(dataSim = dataSim,modelLmer = mod,modelResults = mod.broom ))}
We call the function once and extract the results of this single simulation:
Show the code
out <-FUN_sim_2wayInt(n_subj =1000, # number of subjectsn_subj_prop_p =c(.5, .5), # proportion of republican and democrat subjectsn_subj_prop_e =c(.5, .5), # proportion of subjects without and with extreme weather exposuren_trial =25, # number of trialsbeta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # main effect of political affiliation (polAff)beta_e =0, # main effect of extreme weather exposure (ewe)beta_p_e_inx =0.4/3.5, # two-way interaction effect of polAff and ewesubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible numbers be truncuated?setSeed =1234# seed number to achieve reproducible results. Set to NULL for simulations!)# Get results tableresultsTable <- out$modelResults %>%select(-c(std.error, statistic, df)) %>%mutate(across(where(is_double), ~round(.x, 4))) %>% knitr::kable()formulaUsedForFit <-paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse =" ")# Create plotp.demo.2wayInt <- out$dataSim %>%ggplot(aes(x = ewe, y = deltaDuration, color = polAff)) +geom_hline(yintercept =0) +geom_violin(alpha =0.3) +geom_point(data = polAff_ewe_trueEffects,mapping =aes(x = ewe, y = trueDeltaDuration, fill = polAff),shape ="circle open",size =3.5,stroke =2,color ="black",position =position_dodge(width = .9),show.legend =FALSE ) +stat_summary(fun = mean,fun.min = \(x){mean(x) -sd(x)},fun.max = \(x){mean(x) +sd(x)},position =position_dodge(width = .9) ) + ggrepel::geom_label_repel(data = polAff_ewe_trueEffects,mapping =aes(x = ewe, y = trueDeltaDuration, fill = polAff, label =round(trueDeltaDuration, 4)),color ="black",box.padding =1,position =position_dodge(width = .9),show.legend =FALSE ) +scale_color_manual(values =c("red", "dodgerblue")) +scale_fill_manual(values =c("white", "white")) +scale_y_continuous(breaks =seq(-1, 1, .2)) +labs(title ="Demo Output of One Simulation for Two-Way Interaction") +theme_bw()jpeg(file ="../images/pa_demo2wayInt.jpeg",width =9.1, height =6.3, units ="in", res =600)print(p.demo.2wayInt)invisible(dev.off())
Figure 10 visualizes the results of this single simulation and Table 4 summarizes the statistical results of fitting the actual model used in data generation to the simulated data.
Figure 10: Visual representation of results of one simulation created using FUN_sim_2wayInt. Violin plots display the full distribution of the data. Points and surrounding lines indicate the mean ± 1 SD. The black horizontal line displays the true sample mean and the black open circles indicate the true means for each cell.
Table 4: Statistical results of one simulation created using FUN_sim_2wayInt. Data was fit using deltaDuration ~ polAff * ewe + (1 | subj) + (1 | trial).
effect
group
term
estimate
p.value
fixed
NA
(Intercept)
-0.0006
0.9625
fixed
NA
polAff.dem-rep
0.1290
0.0000
fixed
NA
ewe.TRUE-FALSE
-0.0139
0.4516
fixed
NA
polAff.dem-rep:ewe.TRUE-FALSE
0.1398
0.0002
ran_pars
subj
sd__(Intercept)
0.2855
NA
ran_pars
trial
sd__(Intercept)
0.0382
NA
ran_pars
Residual
sd__Observation
0.3232
NA
3.5.3 Power Simulation
In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run.
Show the code
FUN_sim_2wayInt_pwr <-function(sim, ...){ out <-FUN_sim_2wayInt(...) modelResults <- out$modelResults %>%mutate(sim = sim) %>%relocate(sim)return(modelResults)}# How many simulations should be run?n_sims <-1000# What are the breaks for number of subjects we would like to calculate power for?breaks_subj <-c(900, 950, 1000)# What are the breaks for SESOI?breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)# What are the breaks for different error SDs?breaks_sigma <-c((.29+.04), 2*(.29+.04))res_2wayInt <-tibble()for (s inseq_along(breaks_sigma)) { res_sesoi <-tibble()for (sesoi inseq_along(breaks_sesoi)) { res_nSubj <-tibble()for (nSubj inseq_along(breaks_subj)) {# Give feedback regarding which model is simulatedcat(paste0("Simulation:\n"," sigma = ", round(breaks_sigma[s], 4), "\n"," sesoi = ", round(breaks_sesoi[sesoi], 4), "\n"," nSubject = ", breaks_subj[nSubj], "\n" ))# Start timercat(paste0("Start date time: ", lubridate::now(), "\n"))tic()# Loop over simulations pwr <-map_df(1:n_sims, FUN_sim_2wayInt_pwr,n_subj = breaks_subj[nSubj],beta_p_e_inx = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Stop timer and calculate elapsed time elapsed_time <-toc(quiet =TRUE) elapsed_seconds <- elapsed_time$toc - elapsed_time$tic elapsed_minutes <- elapsed_seconds /60cat(paste0("End date time: ", lubridate::now(), "\n"))cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")# Add number of subjects to pwr pwr <- pwr %>%mutate(nSubjects = breaks_subj[nSubj],sesoi = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Add results to the results table res_nSubj <- res_nSubj %>%rbind(pwr) }# Add results to the results table res_sesoi <- res_sesoi %>%rbind(res_nSubj) }# Add results to the results table res_2wayInt <- res_2wayInt %>%rbind(res_sesoi)}res_2wayInt.summary <- res_2wayInt %>%filter(term =="polAff.dem-rep:ewe.TRUE-FALSE") %>%group_by(sigma, sesoi, nSubjects) %>%summarise(power =mean(p.value <0.05),ci.lower =binom.confint(power*n_sims, n_sims, methods ="exact")$lower,ci.upper =binom.confint(power*n_sims, n_sims, methods ="exact")$upper,.groups ='drop' ) %>%mutate(sigma_fact =factor(format(round(sigma, 4), nsmall =4)),sigma_level =match(sigma_fact, levels(sigma_fact)),sesoi_fact =factor(format(round(sesoi, 4), nsmall =4)),sesoi_level =match(sesoi_fact, levels(sesoi_fact)) )# Save results in a list objecttime <-format(Sys.time(), "%Y%m%d_%H%M")fileName <-paste0("res_2wayInt", "_", time, ".RDS")saveRDS(list(res_2wayInt = res_2wayInt,res_2wayInt.summary = res_2wayInt.summary ),file =file.path("../powerSimulationsOutput", fileName))
We retrieve pre-run results:
Show the code
# Load power simulation dataresList_2wayInt <-readRDS(file.path("../powerSimulationsOutput", "res_2wayInt_20240814_1052.RDS"))resList_2wayInt.summary <- resList_2wayInt$res_2wayInt.summary# Extract power values for some specific assumptionschosenN <-950chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1429", "0.1714")powerValues <- resList_2wayInt.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN) %>%mutate(power_str =paste0(round(power*100, 2), "%")) %>%pull(power_str)# Extract number of simulationslabel_nSimulations <- resList_2wayInt$res_2wayInt$sim %>%n_distinct()# Repeat breaks_sesoibreaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)
Figure 11 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.
Figure 11: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff * ewe + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.
Figure 12 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.
Figure 12: Power curves for the two-way interaction polAff × ewe. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes, taking into account participant exclusion-rates of 10% (N = 900), 5% (N = 950), and 0% (N = 1000). Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.
Show the code
# Get power values of interestchosenN <-950chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1429", "0.1714")powerValues <- resList_2wayInt.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN)# Get lower CI for power values for more liberal and more conservative sigma assumptionspowerValues_sigmaLib <- powerValues %>%filter(sigma_fact =="0.3300")powerValues_sigmaCons <- powerValues %>%filter(sigma_fact =="0.6600")# Interpolate the effect sizes at which we achieve 95% powereff_sigmaLib <-approx(powerValues_sigmaLib$ci.lower, powerValues_sigmaLib$sesoi, xout =0.95)$yeff_sigmaCons <-approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$sesoi, xout =0.95)$y# Round results for display in texteff_sigmaLib_txt <-format(round(eff_sigmaLib, 4), nsmall =4)eff_sigmaCons_txt <-format(round(eff_sigmaCons, 4), nsmall =4)
To assess the smallest effect size that can be detected with 95% statistical power, we inspect the lower bounds of the 95%-CI power estimates in Figure 12. Specifically, we focus on the power simulation results for N = 950, which takes into account a participant exclusion-rate of 5%. There, we interpolate between the two point estimates that lie just below and above the 95% power line, i.e., between the power estimates for effect sizes 0.1429 and 0.1714. Assuming an error SD of 0.3300, we achieve 95% statistical power to detect a two-way interaction effect of at least 0.1530. For a more conservative error SD of 0.6600, this smallest detectable effect size is only marginally higher (0.1619).
3.5.3.1 Increased Sample Size
We additionally conduct an effect-size sensitivity analysis for increased sample sizes.
Show the code
FUN_sim_2wayInt_pwr <-function(sim, ...){ out <-FUN_sim_2wayInt(...) modelResults <- out$modelResults %>%mutate(sim = sim) %>%relocate(sim)return(modelResults)}# How many simulations should be run?n_sims <-1000# What are the breaks for number of subjects we would like to calculate power for?breaks_subj <-c(1500, 1750, 2000)# What are the breaks for SESOI?breaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)# What are the breaks for different error SDs?breaks_sigma <-c((.29+.04), 2*(.29+.04))res_2wayInt <-tibble()for (s inseq_along(breaks_sigma)) { res_sesoi <-tibble()for (sesoi inseq_along(breaks_sesoi)) { res_nSubj <-tibble()for (nSubj inseq_along(breaks_subj)) {# Give feedback regarding which model is simulatedcat(paste0("Simulation:\n"," sigma = ", round(breaks_sigma[s], 4), "\n"," sesoi = ", round(breaks_sesoi[sesoi], 4), "\n"," nSubject = ", breaks_subj[nSubj], "\n" ))# Start timercat(paste0("Start date time: ", lubridate::now(), "\n"))tic()# Loop over simulations pwr <-map_df(1:n_sims, FUN_sim_2wayInt_pwr,n_subj = breaks_subj[nSubj],beta_p_e_inx = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Stop timer and calculate elapsed time elapsed_time <-toc(quiet =TRUE) elapsed_seconds <- elapsed_time$toc - elapsed_time$tic elapsed_minutes <- elapsed_seconds /60cat(paste0("End date time: ", lubridate::now(), "\n"))cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")# Add number of subjects to pwr pwr <- pwr %>%mutate(nSubjects = breaks_subj[nSubj],sesoi = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Add results to the results table res_nSubj <- res_nSubj %>%rbind(pwr) }# Add results to the results table res_sesoi <- res_sesoi %>%rbind(res_nSubj) }# Add results to the results table res_2wayInt <- res_2wayInt %>%rbind(res_sesoi)}res_2wayInt.summary <- res_2wayInt %>%filter(term =="polAffX_p:eweX_e") %>%group_by(sigma, sesoi, nSubjects) %>%summarise(power =mean(p.value <0.05),ci.lower =binom.confint(power*n_sims, n_sims, methods ="exact")$lower,ci.upper =binom.confint(power*n_sims, n_sims, methods ="exact")$upper,.groups ='drop' ) %>%mutate(sigma_fact =factor(format(round(sigma, 4), nsmall =4)),sigma_level =match(sigma_fact, levels(sigma_fact)),sesoi_fact =factor(format(round(sesoi, 4), nsmall =4)),sesoi_level =match(sesoi_fact, levels(sesoi_fact)) )# Save results in a list objecttime <-format(Sys.time(), "%Y%m%d_%H%M")fileName <-paste0("res_2wayInt_increasedN", "_", time, ".RDS")saveRDS(list(res_2wayInt = res_2wayInt,res_2wayInt.summary = res_2wayInt.summary ),file =file.path("../powerSimulationsOutput", fileName))
Show the code
# Load power simulation dataresList_2wayInt <-readRDS(file.path("../powerSimulationsOutput", "res_2wayInt_increasedN_20250308_0703.RDS"))resList_2wayInt.summary <- resList_2wayInt$res_2wayInt.summary# Extract power values for some specific assumptionschosenN <-2000chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1429", "0.1714")powerValues <- resList_2wayInt.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN) %>%mutate(power_str =paste0(round(power*100, 2), "%")) %>%pull(power_str)# Extract number of simulationslabel_nSimulations <- resList_2wayInt$res_2wayInt$sim %>%n_distinct()# Repeat breaks_sesoibreaks_sesoi <- (0.4/3.5)*seq(1, 2, .25)
Figure 12 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.
Figure 13: Power curves for the two-way interaction polAff × ewe. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes. Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.
3.6 Three-Way Interaction Effect
As argued in the Registered Report, we hypothesize:
(H3): The two-way interaction effect of extreme weather exposure and political affiliation is greater for individuals who more strongly attribute extreme weather events to climate change. In other words, there is a positive three-way interaction effect of political affiliation, exposure to extreme weather events, and their attribution to climate change on ΔDuration.
Show the code
# Smallest Effect Size Of Interest (SESOI)SESOI <-0.4/3.5# Betasbeta_p <- SESOIbeta_e <-0beta_s <-0beta_p_e_inx <- SESOIbeta_p_s_inx <-0beta_e_s_inx <-0beta_p_e_s_inx <- SESOI# Predicted "true" effectspolAff_ewe_subjAttr_trueEffects <-expand_grid(polAff =factor(c("rep", "dem"), levels =c("rep", "dem")),ewe =factor(c(FALSE, TRUE), levels =c(FALSE, TRUE)),subjAttr =factor(c("-SD", "M", "+SD"), levels =c("-SD", "M", "+SD"))) %>%add_contrast("polAff", contrast ="anova", colnames ="X_p") %>%add_contrast("ewe", contrast ="anova", colnames ="X_e") %>%add_contrast("subjAttr", contrast ="anova", colnames =c("X_s_1", "X_s_2")) %>%mutate(trueDeltaDuration =0+# intercept X_p * beta_p +# main effect polAff X_e * beta_e +# main effect ewe X_s_1 * beta_s +# main effect subjAttr, dummy variable for M vs. -SD X_s_2 * (2* beta_s) +# main effect subjAttr, dummy variable for +SD vs. -SD X_p * X_e * beta_p_e_inx +# 2-way interaction polAff * ewe X_p * X_s_1 * beta_p_s_inx +# 2-way interaction polAff * subjAttr, dummy variable for M vs. -SD X_p * X_s_2 * (2* beta_p_s_inx) +# 2-way interaction polAff * subjAttr, dummy variable for +SD vs. -SD X_e * X_s_1 * beta_e_s_inx +# 2-way interaction ewe * subjAttr, dummy variable for M vs. -SD X_e * X_s_2 * (2* beta_e_s_inx) +# 2-way interaction ewe * subjAttr, dummy variable for +SD vs. -SD X_p * X_e * X_s_1 * beta_p_e_s_inx +# 3-way interaction polAff*ewe*subjAttr, dummy variable for M vs -SD X_p * X_e * X_s_2 * (2* beta_p_e_s_inx) # 3-way interaction polAff*ewe*subjAttr, dummy variable for +SD vs. -SD )
3.6.1 SESOI for Three-way Interaction
In the sections above, we derived the SESOI used for the main effect sample-size determination analysis and for the two-way interaction effect-size sensitivity analysis for the binary variables polAff (dem vs. rep) and ewe (TRUE vs. FALSE). subjAttr, however, is a continuous variable that ranges from 1 to 5. Fortunately, in finding a theoretically sound SESOI for the three-way interaction, the same considerations apply as for the main effect and two-way interaction before. We just need to translate these considerations into the continuous metric of subjAttr.
We start by noticing that the complete fixed three-way interaction polAff × ewe × subjAttr is modeled as:
Now, let’s calculate this two-way interaction for two individuals who differ in their level of subjective attribution of extreme weather events to climate change. First, an individual who has an average score on subjAttr will show the following two-way interaction effect, with \(\mu_{subjAttr}\) being the sample average of the variable subjAttr:
Second, we define an individual with a low score on subjAttr as one that shows a subjective attribution of one SD bellow the average. This individual will show the following two-way interaction effect, with \(\sigma_{subjAttr}\) being the SD of subjAttr:
As outlined above in Section Section 3.5.1, we assume that the SESOI for the two-way interaction polAff × ewe is 0.1143 for an average individual (with respect to subjAttr). If the same two-way interaction polAff × ewe shrinks to zero for an individual low in subjAttr, we would consider this interaction effect difference as theoretically relevant (see Figure 14). These assumptions translate to:
We will assess subjective attribution of extreme weather events to climate change using the same questions, response options, and aggregation as Ogunbode et al. (2019). These authors reported \(\mu_{subjAttr}\) = 3.67 and \(\sigma_{subjAttr}\) = 0.85, resulting in:
Figure 14: Assumed ΔDuration values for individuals scoring low (-SD) and average (M) on subjAttr. While individuals scoring average on subjAttr show a two-way interaction effect of polAff × ewe of 0.1143, this effect shrinks to zero for individuals scoring low on subjAttr. These individuals only show the predicted main effect of polAff (whis is 0.1143).
3.6.2 Data Simulation Function
We finally define a function that simulates data for the three-way interaction effect of political affiliation with extreme weather exposure and attribution of extreme weather events to climate change on ΔDuration: FUN_sim_3wayInt. The function will simulate data according to the following model, expressed in lme4-lingo:
The function FUN_sim_3wayInt takes, among others, the following important arguments (in addition to the arguments discussed for FUN_sim_2wayInt):
s_mean: Sample mean of subjective attribution. Based on results reported by Ogunbode et al. (2019), we set this value to 3.67.
s_sd: Sample SD of subjective attribution. Based on results reported by Ogunbode et al. (2019), we set this value to 0.85.
beta_s: Fixed main effect of subjective attribution. As we are interested in the moderating role of subjective attribution of extreme weather events to climate change, we set this main effect to zero.
beta_p_s_inxFixed two-way interaction effect of political affiliation and subjective attribution. In order to accurately model a three-way interaction, one needs to include all two-way interactions in statistical models. Therefore, we include this two-way interaction, but we assume it to be zero.
beta_e_s_inxFixed two-way interaction effect of extreme weather exposure and subjective attribution. As before, we include this interaction for accurately modelling the three-way interaction of interest, but we assume this two-way interaction to be zero.
beta_p_e_s_inxFixed three-way interaction effect of political affiliation, extreme weather exposure, and subjective attribution. We set this initial value to the SESOI derived above, i.e. 0.1345. We investigate how changing this effect size impacts statistical power, as we are conducting effect-size sensitivity analyses for interaction effects.
The function is defined below:
Show the code
# define data simulation functionFUN_sim_3wayInt <-function(n_subj =1000, # number of subjectsn_subj_prop_p =c(.5, .5), # proportion of republican and democrat subjectsn_subj_prop_e =c(.5, .5), # proportion of subjects without and with extreme weather exposuren_subj_prop_s =c(.5, .5), # proportion of subjects with low and high subjective attribution of EWE to CCn_trial =25, # number of trialss_mean =3.67, # mean of subjAttr, see Ogunbode et al. (2019)s_sd =0.85, # sd of subjAttr, see Ogunbode et al. (2019)beta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # main effect of political affiliation (polAff)beta_e =0, # main effect of extreme weather exposure (ewe)beta_s =0, # main effect of subjective attribution of ewe to climate change (subjAttr)beta_p_e_inx =0.4/3.5, # two-way interaction effect of polAff and ewebeta_p_s_inx =0, # two-way interaction effect of polAff and subjAttrbeta_e_s_inx =0, # two-way interaction effect of ewe and subjAttrbeta_p_e_s_inx = (0.4/3.5)/0.85, # three-way interaction effect of polAff, ewe, and subjAttrsubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible numbers be truncuated?setSeed =NULL# seed number to achieve reproducible results. Set to NULL for simulations!) {# set seed to achieve reproducible results for demonstration purposesset.seed(setSeed)# simulate data for dwell time on carbon information dataSim <-# add random factor subjectadd_random(subj = n_subj) %>%# add random factor trialadd_random(trial = n_trial) %>%# add between-subject factor political affiliation (with anova contrast)add_between("subj", polAff =c("rep", "dem"), .prob = n_subj_prop_p*n_subj, .shuffle =TRUE) %>%add_contrast("polAff", colnames ="X_p", contrast ="anova") %>%# add between-subject factor extreme weather exposure (with anova contrast)add_between("subj", ewe =c(FALSE, TRUE), .prob = n_subj_prop_e*n_subj, .shuffle =TRUE) %>%add_contrast("ewe", colnames ="X_e", contrast ="anova") %>%# add between-subject variable subjective attribution of EWE to climate changemutate(subjAttr =rep(rnorm(n = n_subj, mean = s_mean, sd = s_sd), each = n_trial),subjAttr_c =scale(subjAttr, center =TRUE, scale =FALSE)[,1] ) %>%# add by-subject random interceptadd_ranef("subj", S_0 = subj_0) %>%# add by-trial random interceptadd_ranef("trial", T_0 = trial_0) %>%# add error termadd_ranef(e_st = sigma) %>%# add response valuesmutate(# add together fixed and random effects for each effectB_0 = beta_0 + S_0 + T_0,B_p = beta_p,B_e = beta_e,B_s = beta_s,B_p_e_inx = beta_p_e_inx,B_p_s_inx = beta_p_s_inx,B_e_s_inx = beta_e_s_inx,B_p_e_s_inx = beta_p_e_s_inx,# calculate dv by adding each effect term multiplied by the relevant# effect-coded factors and adding the error termdeltaDuration = B_0 + e_st + (X_p * B_p) + (X_e * B_e) + (subjAttr_c * B_s) + (X_p * X_e * B_p_e_inx) + (X_p * subjAttr_c * B_p_s_inx) + (X_e * subjAttr_c * B_e_s_inx) + (X_p * X_e * subjAttr_c * B_p_e_s_inx) )# unset seedset.seed(NULL)# truncuate impossible deltaDurationsif(truncNums) { dataSim <- dataSim %>%mutate(deltaDuration =if_else(deltaDuration <-1, -1,if_else(deltaDuration >1, 1, deltaDuration))) }# run a linear mixed effects model and check summary mod <-lmer( deltaDuration ~ polAff*ewe*subjAttr_c + (1| subj) + (1| trial),data = dataSim ) mod.sum <-summary(mod)# get results in tidy format mod.broom <- broom.mixed::tidy(mod)return(list(dataSim = dataSim,modelLmer = mod,modelResults = mod.broom ))}
We call the function once and extract the results of this single simulation:
Show the code
# Note to myself: Consider setting beta_p = 0.4/3.5, beta_p_e_inx = 0.4/3.5,# and beta_p_e_s_inx to 2 * 0.4/3.5/(2*.85).# This way, the interaction effect polAff:ewe is 0.4/3.5 for individuals# with an average subjAttr (subjAttr_c = 0). For individuals with# subjAttr = mean - SD, polAff:ewe is 0. For individuals with subjAttr = mean + SD,# polAff:ewe is 2 * 0.4/3.5.out <-FUN_sim_3wayInt(n_subj =1000, # number of subjectsn_subj_prop_p =c(.5, .5), # proportion of republican and democrat subjectsn_subj_prop_e =c(.5, .5), # proportion of subjects without and with extreme weather exposuren_subj_prop_s =c(.5, .5), # proportion of subjects with low and high subjective attribution of EWE to CCn_trial =25, # number of trialss_mean =3.67, # mean of subjAttr, see Ogunbode et al. (2019)s_sd =0.85, # sd of subjAttr, see Ogunbode et al. (2019)beta_0 =0, # intercept (grand mean) for deltaDurationbeta_p =0.4/3.5, # main effect of political affiliation (polAff)beta_e =0, # main effect of extreme weather exposure (ewe)beta_s =0, # main effect of subjective attribution of ewe to climate change (subjAttr)beta_p_e_inx =0.4/3.5, # two-way interaction effect of polAff and ewebeta_p_s_inx =0, # two-way interaction effect of polAff and subjAttrbeta_e_s_inx =0, # two-way interaction effect of ewe and subjAttrbeta_p_e_s_inx = (0.4/3.5)/0.85, # three-way interaction effect of polAff, ewe, and subjAttrsubj_0 = .29, # by-subject random intercept sd for dt carbontrial_0 = .04, # by-trial random intercept sdsigma =1*(.29+.04), # residual (error) sdtruncNums =TRUE, # should impossible numbers be truncuated?setSeed =123# seed number to achieve reproducible results. Set to NULL for simulations!)# Get results tableresultsTable <- out$modelResults %>%select(-c(std.error, statistic, df)) %>%mutate(across(where(is_double), ~round(.x, 4))) %>% knitr::kable()formulaUsedForFit <-paste(as.character(formula(out$modelLmer))[c(2,1,3)], collapse =" ")# Create predictions plot# refit model to dispaly subjAttr levels in original metric (not mean centered)m <-lmer( deltaDuration ~ polAff*ewe*subjAttr + (1| subj) + (1| trial),data = out$dataSim)# define spotlights for spotlight analysisspotlights <-c(3.67-0.85, 3.67, 3.67+0.85)# create plot showing predictionsp.demo.3wayInt.pred <-predict_response(m, terms =c("ewe", "polAff", "subjAttr[spotlights]"))p.demo.3wayInt <- p.demo.3wayInt.pred %>%plot(colors =c("red", "dodgerblue")) +coord_cartesian(ylim =c(-.25, .25)) +theme_bw()jpeg(file ="../images/pa_demo3wayInt.jpeg",width =9.1, height =6.3, units ="in", res =600)print(p.demo.3wayInt)invisible(dev.off())
Figure 15 visualizes predictions based on this single simulation and Table 5 summarizes the statistical results of fitting the actual model used in data generation to the simulated data.
Figure 15: Visual representation of results of one simulation created using FUN_sim_3wayInt. Points indicate the predicted means surrounded by 95% Confidence Intervals. Panels indicate predictions for different values of subjAttr (Mean - SD, Mean, and Mean + SD).
Table 5: Statistical results of one simulation created using FUN_sim_3wayInt. Data was fit using deltaDuration ~ polAff * ewe * subjAttr_c + (1 | subj) + (1 | trial). Note that subjAttr is mean centered to ease interpretation of lower-level interactions and main effects.
effect
group
term
estimate
p.value
fixed
NA
(Intercept)
-0.0022
0.8478
fixed
NA
polAff.dem-rep
0.1149
0.0000
fixed
NA
ewe.TRUE-FALSE
-0.0023
0.9031
fixed
NA
subjAttr_c
0.0091
0.4141
fixed
NA
polAff.dem-rep:ewe.TRUE-FALSE
0.1426
0.0001
fixed
NA
polAff.dem-rep:subjAttr_c
0.0187
0.4025
fixed
NA
ewe.TRUE-FALSE:subjAttr_c
0.0142
0.5239
fixed
NA
polAff.dem-rep:ewe.TRUE-FALSE:subjAttr_c
0.1111
0.0130
ran_pars
subj
sd__(Intercept)
0.2849
NA
ran_pars
trial
sd__(Intercept)
0.0323
NA
ran_pars
Residual
sd__Observation
0.3240
NA
3.6.3 Power Simulation
In the following code, the simulations are calculated. We do not recommend executing this code junk as it takes several hours to run.
Show the code
FUN_sim_3wayInt_pwr <-function(sim, ...){ out <-FUN_sim_3wayInt(...) modelResults <- out$modelResults %>%mutate(sim = sim) %>%relocate(sim)return(modelResults)}# How many simulations should be run?n_sims <-1000# What are the breaks for number of subjects we would like to calculate power for?breaks_subj <-c(900, 950, 1000)# What are the breaks for SESOI?breaks_sesoi <- (0.4/3.5)/0.85*seq(1, 2, .25)# What are the breaks for different error SDs?breaks_sigma <-c((.29+.04), 2*(.29+.04))res_3wayInt <-tibble()for (s inseq_along(breaks_sigma)) { res_sesoi <-tibble()for (sesoi inseq_along(breaks_sesoi)) { res_nSubj <-tibble()for (nSubj inseq_along(breaks_subj)) {# Give feedback regarding which model is simulatedcat(paste0("Simulation:\n"," sigma = ", round(breaks_sigma[s], 4), "\n"," sesoi = ", round(breaks_sesoi[sesoi], 4), "\n"," nSubject = ", breaks_subj[nSubj], "\n" ))# Start timercat(paste0("Start date time: ", lubridate::now(), "\n"))tic()# Loop over simulations pwr <-map_df(1:n_sims, FUN_sim_3wayInt_pwr,n_subj = breaks_subj[nSubj],beta_p_e_s_inx = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Stop timer and calculate elapsed time elapsed_time <-toc(quiet =TRUE) elapsed_seconds <- elapsed_time$toc - elapsed_time$tic elapsed_minutes <- elapsed_seconds /60cat(paste0("End date time: ", lubridate::now(), "\n"))cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")# Add number of subjects to pwr pwr <- pwr %>%mutate(nSubjects = breaks_subj[nSubj],sesoi = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Add results to the results table res_nSubj <- res_nSubj %>%rbind(pwr) }# Add results to the results table res_sesoi <- res_sesoi %>%rbind(res_nSubj) }# Add results to the results table res_3wayInt <- res_3wayInt %>%rbind(res_sesoi)}res_3wayInt.summary <- res_3wayInt %>%filter(term =="polAff.dem-rep:ewe.TRUE-FALSE:subjAttr_c") %>%group_by(sigma, sesoi, nSubjects) %>%summarise(power =mean(p.value <0.05),ci.lower =binom.confint(power*n_sims, n_sims, methods ="exact")$lower,ci.upper =binom.confint(power*n_sims, n_sims, methods ="exact")$upper,.groups ='drop' ) %>%mutate(sigma_fact =factor(format(round(sigma, 4), nsmall =4)),sigma_level =match(sigma_fact, levels(sigma_fact)),sesoi_fact =factor(format(round(sesoi, 4), nsmall =4)),sesoi_level =match(sesoi_fact, levels(sesoi_fact)) )# Save results in a list objecttime <-format(Sys.time(), "%Y%m%d_%H%M")fileName <-paste0("res_3wayInt", "_", time, ".RDS")saveRDS(list(res_3wayInt = res_3wayInt,res_3wayInt.summary = res_3wayInt.summary ),file =file.path("../powerSimulationsOutput", fileName))
We retrieve pre-run results:
Show the code
# Load power simulation dataresList_3wayInt <-readRDS(file.path("../powerSimulationsOutput", "res_3wayInt_20240814_1612.RDS")) resList_3wayInt.summary <- resList_3wayInt$res_3wayInt.summary# Extract power values for some specific effect sizes at N = 1000powerValues <- resList_3wayInt.summary %>%filter(sigma_fact =="0.3300") %>%filter(sesoi_fact =="0.1345") %>%filter(nSubjects ==950) %>%mutate(power_str =paste0(round(power*100, 2), "%")) %>%pull(power_str)# Extract number of simulationslabel_nSimulations <- resList_3wayInt$res_3wayInt$sim %>%n_distinct()# Repeat breaks_sesoibreaks_sesoi <- (0.4/3.5)/(0.85) *seq(1, 2, .25)
Figure 16 displays the distribution of estimated fixed effects across all simulations. The figure shows that the estimated fixed effects are close to the true ones provided as input in the data simulation function, validating that simulations worked as expected.
Figure 16: Distribution of estimated fixed effects resulting from 1000 simulations for the model deltaDuration ~ polAff * ewe * subjAttr_c + (1 | subj) + (1 | trial). Shaded area represent densities, annotated points indicate medians, and thick and thin lines represent 66% and 95% quantiles.
Figure 17 shows results of our effect-size sensitivity analyses. We plot statistical power (y-axis) for different effect sizes (x-axis), taking into account different assumptions for the error SD (color) and sample size (panel). Regarding the latter, we report results not only for the full sample size we aim for (N = 1000), but also for sample sizes taking into account different participant exclusion-rates due to exclusion criteria defined in the Registered Report.
Figure 17: Power curves for the three-way interaction polAff × ewe × subjAttr. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes, taking into account participant exclusion-rates of 10% (N = 900), 5% (N = 950), and 0% (N = 1000). Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.
Show the code
# Get power values of interestchosenN <-950chosenSigma <-c("0.3300", "0.6600")chosenSESOI <-c("0.1681", "0.2017")powerValues <- resList_3wayInt.summary %>%filter(sigma_fact %in% chosenSigma) %>%filter(sesoi_fact %in% chosenSESOI) %>%filter(nSubjects %in% chosenN)# Get lower CI for power values for more liberal and more conservative sigma assumptionspowerValues_sigmaLib <- powerValues %>%filter(sigma_fact =="0.3300")powerValues_sigmaCons <- powerValues %>%filter(sigma_fact =="0.6600")# Interpolate the effect sizes at which we achieve 95% powereff_sigmaLib <-approx(powerValues_sigmaLib$ci.lower, powerValues_sigmaLib$sesoi, xout =0.95)$yeff_sigmaCons <-approx(powerValues_sigmaCons$ci.lower, powerValues_sigmaCons$sesoi, xout =0.95)$y# Round results for display in texteff_sigmaLib_txt_3way <-format(round(eff_sigmaLib, 4), nsmall =4)eff_sigmaCons_txt_3way <-format(round(eff_sigmaCons, 4), nsmall =4)
To assess the smallest effect size that can be detected with 95% statistical power, we inspect the lower bounds of the 95%-CI power estimates in Figure 17. Specifically, we focus on the power simulation results for N = 950, which takes into account a participant exclusion rate of 5%. There, we interpolate between the two point estimates that lie just below and above the 95% power line, i.e., between the power estimates for effect sizes 0.1681 and 0.2017. Assuming an error SD of 0.3300, we achieve 95% statistical power to detect a three-way interaction effect of at least 0.1728. For a more conservative error SD of 0.6600, this smallest detectable effect size is only marginally higher (0.1890).
3.6.3.1 Increased Sample Size
We additionally conduct an effect-size sensitivity analysis for increased sample sizes.
Show the code
FUN_sim_3wayInt_pwr <-function(sim, ...){ out <-FUN_sim_3wayInt(...) modelResults <- out$modelResults %>%mutate(sim = sim) %>%relocate(sim)return(modelResults)}# How many simulations should be run?n_sims <-1000# What are the breaks for number of subjects we would like to calculate power for?breaks_subj <-c(1500, 1750, 2000)# What are the breaks for SESOI?breaks_sesoi <- (0.4/3.5)/0.85*seq(1, 2, .25)# What are the breaks for different error SDs?breaks_sigma <-c((.29+.04), 2*(.29+.04))res_3wayInt <-tibble()for (s inseq_along(breaks_sigma)) { res_sesoi <-tibble()for (sesoi inseq_along(breaks_sesoi)) { res_nSubj <-tibble()for (nSubj inseq_along(breaks_subj)) {# Give feedback regarding which model is simulatedcat(paste0("Simulation:\n"," sigma = ", round(breaks_sigma[s], 4), "\n"," sesoi = ", round(breaks_sesoi[sesoi], 4), "\n"," nSubject = ", breaks_subj[nSubj], "\n" ))# Start timercat(paste0("Start date time: ", lubridate::now(), "\n"))tic()# Loop over simulations pwr <-map_df(1:n_sims, FUN_sim_3wayInt_pwr,n_subj = breaks_subj[nSubj],beta_p_e_s_inx = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Stop timer and calculate elapsed time elapsed_time <-toc(quiet =TRUE) elapsed_seconds <- elapsed_time$toc - elapsed_time$tic elapsed_minutes <- elapsed_seconds /60cat(paste0("End date time: ", lubridate::now(), "\n"))cat("Elapsed time: ", elapsed_minutes, " minutes\n\n")# Add number of subjects to pwr pwr <- pwr %>%mutate(nSubjects = breaks_subj[nSubj],sesoi = breaks_sesoi[sesoi],sigma = breaks_sigma[s] )# Add results to the results table res_nSubj <- res_nSubj %>%rbind(pwr) }# Add results to the results table res_sesoi <- res_sesoi %>%rbind(res_nSubj) }# Add results to the results table res_3wayInt <- res_3wayInt %>%rbind(res_sesoi)}res_3wayInt.summary <- res_3wayInt %>%filter(term =="polAffX_p:eweX_e:subjAttr_c") %>%group_by(sigma, sesoi, nSubjects) %>%summarise(power =mean(p.value <0.05),ci.lower =binom.confint(power*n_sims, n_sims, methods ="exact")$lower,ci.upper =binom.confint(power*n_sims, n_sims, methods ="exact")$upper,.groups ='drop' ) %>%mutate(sigma_fact =factor(format(round(sigma, 4), nsmall =4)),sigma_level =match(sigma_fact, levels(sigma_fact)),sesoi_fact =factor(format(round(sesoi, 4), nsmall =4)),sesoi_level =match(sesoi_fact, levels(sesoi_fact)) )# Save results in a list objecttime <-format(Sys.time(), "%Y%m%d_%H%M")fileName <-paste0("res_3wayInt_increasedN", "_", time, ".RDS")saveRDS(list(res_3wayInt = res_3wayInt,res_3wayInt.summary = res_3wayInt.summary ),file =file.path("../powerSimulationsOutput", fileName))
Show the code
# Load power simulation dataresList_3wayInt <-readRDS(file.path("../powerSimulationsOutput", "res_3wayInt_increasedN_20250308_1649.RDS"))resList_3wayInt.summary <- resList_3wayInt$res_3wayInt.summary# Extract power values for some specific effect sizes at N = 1000powerValues <- resList_3wayInt.summary %>%filter(sigma_fact =="0.3300") %>%filter(sesoi_fact =="0.1345") %>%filter(nSubjects ==950) %>%mutate(power_str =paste0(round(power*100, 2), "%")) %>%pull(power_str)# Extract number of simulationslabel_nSimulations <- resList_3wayInt$res_3wayInt$sim %>%n_distinct()# Repeat breaks_sesoibreaks_sesoi <- (0.4/3.5)/(0.85) *seq(1, 2, .25)
Figure 18: Power curves for the three-way interaction polAff × ewe × subjAttr. Points represent simulated power surrounded by a 95%-CI based on 1000 simulations with α = 0.05. Note that, in contrast to Figure 9, the x-axsis represents different effect sizes, starting from the defined SESOI, while the panels represent different sample sizes. Note that estimates are displayed with a slight shift along the x-axis to reduce overlap. Since we only consider effect sizes equal to or greater than the SESOI to be of relevance, power curves were only examined for positive effect sizes, which is why only positive values are shown.
3.7 Conclusion
Using power simulations with mixed-effects models for sample-size determination and effect-size sensitivity analyses, we show that with a final sample size of N = 950:
We will achieve \(\ge\) 95% statistical power to detect a main effect of political affiliation.
We will be able to detect a two-way interaction effect of political affiliation with extreme weather exposure of at least 0.1619 with \(\ge\) 95% statistical power.
We will be able to detect a three-way interaction effect of political affiliation with extreme weather exposure and subjective attribution of extreme weather events to climate change of at least 0.1890 with \(\ge\) 95% statistical power.
By increasing the planned sample size to N = 2’000, we will be able to detect all of the defined smallest effect sizes of interest of the hypothesized main effect and the two- and three-way interaction effects with a statistical power of at least 95%.
Berger, Sebastian, and Annika M Wyss. 2021. “Measuring Pro-Environmental Behavior Using the Carbon Emission Task.”Journal of Environmental Psychology 75: 101613. https://doi.org/https://doi.org/10.1016/j.jenvp.2021.101613.
Giner-Sorolla, Roger, Amanda K. Montoya, Alan Reifman, Tom Carpenter, Neil A. Lewis, Christopher L. Aberson, Dries H. Bostyn, et al. 2024. “Power to Detect What? Considerations for Planning and Evaluating Sample Size.”Personality and Social Psychology Review, February, 10888683241228328. https://doi.org/10.1177/10888683241228328.
IPCC, ed. 2023. “Weather and Climate Extreme Events in a Changing Climate.” In, 1513–1766. Cambridge: Cambridge University Press. https://doi.org/10.1017/9781009157896.013.
Konisky, David M., Llewelyn Hughes, and Charles H. Kaylor. 2016. “Extreme Weather Events and Climate Change Concern.”Climatic Change 134 (4): 533–47. https://doi.org/10.1007/s10584-015-1555-3.
Ogunbode, Charles A., Christina Demski, Stuart B. Capstick, and Robert G. Sposato. 2019. “Attribution Matters: Revisiting the Link Between Extreme Weather Experience and Climate Change Mitigation Responses.”Global Environmental Change 54 (January): 31–39. https://doi.org/10.1016/j.gloenvcha.2018.11.005.
Reeck, Crystal, Daniel Wall, and Eric J. Johnson. 2017. “Search Predicts and Changes Patience in Intertemporal Choice.”Proceedings of the National Academy of Sciences 114 (45): 11890–95. https://doi.org/10.1073/pnas.1707040114.
Willemsen, Martijn C., and Eric J. Johnson. 2019. “(Re)visiting the Decision Factory: Obswerving Cognition with MouselabWEB.” In, edited by Michael Schulte-Mecklenbeck, Anton Kuehberger, and Joseph G. Johnson, 2nd ed., 76–95. New York: Routledge. https://doi.org/10.4324/9781315160559.